dat1 <- c(71.9, 0.2, 4.2, 1.5, 2.0, 1.2, 1.1, 1.1, 0.4, 0.7, 0.3, 1.3, 57.1, 2.5, 6.3, 1.7, 1.2, 13.8,
0.4, 3.1, 3.6, 0.7, 139.5, 31.8, 0.1, 0.3, 5.7, 9.8, 3.8, 1.4, 0.5, 0.1, 0.5, 0.4, 3.1, 5.6, 0.1, 0.1,
1.3, 0.05, 37.3, 0.2, 1.2, 1.8, 0.4, 0.1, 2.6, 3.3, 0.4, 12.1, 35.6, 1.7, 10.1, 1.0, 3.1, 0.1, 14.8,
0.05, 0.3, 2.0, 0.05, 3.6, 0.6, 13.5, 0.6, 1.0, 3.7, 0.3, 0.2, 0.1, 0.05, 15.1, 1.7, 23.2, 3.8, 17.8,
6.0, 0.7, 4.1, 0.1, 0.6, 0.05, 0.3, 0.5, 0.1, 0.3, 31.9, 0.05, 0.1, 13.5, 0.2, 0.1, 0.05, 0.2, 0.05,
4.7, 6.6, 2.6, 0.6, 23.0, 0.6, 0.3, 0.1, 4.0, 0.05, 0.05, 0.5, 24.0, 5.7, 1.3, 0.3, 0.3, 0.2, 18.4, 0.05,
0.2, 1.1, 0.1, 2.8, 0.5, 0.1, 0.2, 0.1, 1.6, 2.2, 0.05, 11.5, 1.7, 0.2, 0.05, 0.05, 7.2, 0.3, 0.05,
0.1, 0.1, 3.8, 2.4, 4.2, 1.1, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.05, 0.2, 0.05, 1.0, 0.2, 17.5, 0.2,
0.1, 3.1, 54.3, 0.05, 0.1, 1.0, 0.2, 0.8, 0.05, 0.05, 1.4, 0.05, 6.3, 0.5, 6.2, 0.2, 0.2, 1.2, 0.1, 0.1,
0.8, 3.6, 0.1, 0.1, 0.2, 6.8, 0.05, 1.2, 0.05, 0.05, 6.4, 0.05, 0.3, 0.05, 0.05, 0.05, 0.4, 0.05, 1.1,
0.05, 0.2, 4.0, 0.05, 16.3, 9.6, 0.1, 15.5, 7.5, 0.3, 2.2, 81.0, 1.6, 10.0, 20.6)

dat2 <- c(0, 10, 10, 20, 20, 20, 30, 30, 30, 30, 30, 40, 40, 50, 50, 60, 60, 70)

dat <- list(dat1, dat2)

#return(-ks.test(data, y="pgamma", par[1], par[2])$p.value)

gammaFit <- function(data, par) {
	if(par[1] > 0 && par[2] > 0) {
		xs <- data$mids
		ys <- data$density
		res <- sum((dgamma(xs, par[1], par[2]) - ys)^2)
		return(res)
	}
	Inf
}

normFit <- function(data, par) {
	xs <- data$mids
	ys <- data$density
	res <- sum((dnorm(xs, par[1], par[2]) - ys)^2)
	return(res)
}

normFit2 <- function(data, par) {
	xs <- data$mids
	ys <- data$density
	
	resn <- optim(par = c(E, V), normFit, data = d)
	
	res <- sum(par[1] (dnorm(xs, par[1], par[2]) - ys)^2)
	return(res)
}

par(mfrow=c(1,2))

for(i in 1:2)
{
	da <- dat[[i]]
	E <<- mean(da)
	V <<- var(da)

	d <- hist(dat[[i]], plot=FALSE)
	resn <- optim(par = c(E, V), normFit, data = d)
	resg <- optim(par = c((E^2)/V, E/V), gammaFit, data = d)

	resn2 <- optim(par = c(0.5, 0.5), normFit, data = d)


	title <- sprintf('N(%.0f,%.0f)=%.4f / G(%.0f,%.0f)=%.4f / LN(%.0f,%.0f)=%.4f', resn$par[1], resn$par[2], resn$value, resg$par[1], resg$par[2], resg$value, resln$par[1], resln$par[2], resln$value)
	hist(dat[[i]], freq=FALSE, main=title)
	
	lines(dnorm(1:100, resn$par[1], resn$par[2]), col='red', lwd=2)
	lines(dgamma(1:100, resg$par[1], resg$par[2]), col='blue', lwd=2)
	if(i == 2)
		legend('topright', c('norm', 'gamma'), lwd=c(2,2), col=c('red','blue'))
}

